Datenanalyse / Statistik mit R leicht gemacht!

Data Science Meetup Darmstadt

Datenbasis

Beschreibung im Fließtext

Hier werden R-Ausdrücke direkt in den Text eingebettet. Anwendungsbeispiel: Kontinuierliche Datenerhebung mit sich ändernden Fallzahlen. Die Zahlen im Text werden bei der Berichtserstellung automatisch aktualisiert.

Die folgenden Daten beruhen auf 3000 Alben von 1002 verschiedenen Bands / Künstlern. Sie stammen von der Webseite https://tsort.info/. Wir verwenden die Version 2-8-0044.

Interaktive Daten-Tabelle

DT::datatable(topalbums)

# Zum Ausprobieren: Argument filter; z. B. "filter = "top" oder filter = "bottom"
# DT::datatable(topalbums, filter = "xxx")

# Zum Ausprobieren: Argument options = list(pageLength = x)
# DT::datatable(topalbums, options = list(pageLength = 7))

Gruppenvergleiche

t-Test: Standard-Ausgabe

Rs Standard-Ausgabe ist nicht direkt für Publikationen verwendbar …

my_bands <- topalbums %>% 
  filter(artist %in% c("Bob Dylan", "The Beatles"))

my_test <- t.test(final_score ~ artist,
   alternative = 'two.sided',
   conf.level = .95, var.equal = FALSE,
   data = my_bands)

my_test

    Welch Two Sample t-test

data:  final_score by artist
t = -2.6999, df = 38.571, p-value = 0.01024
alternative hypothesis: true difference in means between group Bob Dylan and group The Beatles is not equal to 0
95 percent confidence interval:
 -10.54416  -1.51018
sample estimates:
  mean in group Bob Dylan mean in group The Beatles 
                 8.175923                 14.203094 

t-Test: Text-Interpretation mit report

Ergebnisse wie oben textlich fundiert und publikationsfähig zu interpretieren, ist nicht für jedeN der schönste Teil der Datenanalyse. Praktisch, wenn es einfach geht:

# R-Paket: report

report(my_test)

The Welch Two Sample t-test testing the difference of final_score by artist (mean in group Bob Dylan = 8.18, mean in group The Beatles = 14.20) suggests that the effect is negative, statistically significant, and large (difference = -6.03, 95% CI [-10.54, -1.51], t(38.57) = -2.70, p = 0.010; Cohen's d = -0.87, 95% CI [-1.52, -0.20])

t-Test: Grafische Darstellung mit ggstatsplot

# R-Paket: ggstatsplot

ggbetweenstats(my_bands,
  x = artist, y = final_score,
  title = "Gesamtpunkte im Vergleich",
  caption = "Datenquelle: tsort.info")

• Boxplot
• Violin-Plot
• Einzelne Datenpunkte
• Angereichert mit statistischen Kennzahlen

Gruppenvergleich: Tabellarisch mit gtsummary

add_p() wählt automatisch einen Test aus; hier: Wilcoxon Rank-Sum Test.

my_bands %>% 
  select(artist, final_score) %>% 
  tbl_summary(by = artist) %>% 
  add_p()

Characteristic	Bob Dylan, N = 261	The Beatles, N = 321	p-value2
final_score	7 (6, 10)	9 (6, 17)	0.049
1 Median (IQR)
2 Wilcoxon rank sum exact test

# Alternative: Hier direkt die Bands ändern

topalbums %>%
  select(artist, final_score) %>% 
  filter(artist %in% c("Madonna", "Bruce Springsteen")) %>%
  tbl_summary(by = artist) %>%
  add_p()

Characteristic	Bruce Springsteen, N = 211	Madonna, N = 201	p-value2
final_score	9.1 (5.7, 13.0)	10.1 (7.4, 14.9)	0.2
1 Median (IQR)
2 Wilcoxon rank sum exact test

Stattdessen t-Test vorgeben, siehe ?add_p():

topalbums %>% 
  filter(artist %in% c("Madonna", "Bruce Springsteen")) %>%
  select(artist, final_score) %>% 
  tbl_summary(by = artist) %>% 
  add_p(
    # perform t-test for all variables
    test = everything() ~ "t.test",
    # assume equal variance in the t-test
    test.args = all_tests("t.test") ~ list(var.equal = TRUE)
)

Characteristic	Bruce Springsteen, N = 211	Madonna, N = 201	p-value2
final_score	9.1 (5.7, 13.0)	10.1 (7.4, 14.9)	0.5
1 Median (IQR)
2 Two Sample t-test

Gruppenvergleich: Vorschläge für andere Tests

Für Fortgeschrittene zum Üben: Andere Fragestellungen testen.

Vorschlag: Unterscheiden sich Alben, die über 5 Punkte erreicht haben, von Alben, die max. 5 Punkte erreicht haben, im durchschnittlichen Alter? Vermutung / Hypothese: Erfolgreichere Alben sind durchschnittlich älter.

Alternative Idee: Unterscheiden sich Alben, die bis zu dem einschneidenden Jahr 1990 veröffentlicht wurden, in der durchschnittlichen Gesamtpunktzahl von Alben, die nach 1990 veröffentlicht wurden?

Oder andere Bands auswählen.

Oder mit anderen Daten arbeiten: topsongs bzw. für aktuellere Musik albums2000, songs2000.

Wer mag, kann, analog zu oben, Test-Einstellungen variieren (t-Test, gleiche Varianzen annehmen ja oder nein).

# Unterscheiden sich Alben, die über 5 Punkte erreicht haben, von Alben, die max. 5 Punkte erreicht haben, im durchschnittlichen Alter?
# Vermutung / Hypothese: Erfolgreichere Alben sind durchschnittlich älter.

topalbums <- topalbums %>% 
  mutate(ueber5 = ifelse(final_score > 5, "> 5 P", "< 5 P"))

# Variable in erster Zeile einfügen
# my_test <- t.test(year ~ xxx,
#    alternative = 'greater',
#    conf.level = .95, var.equal = FALSE,
#    data = topalbums)
# 
# my_test

# Text-Interpretation: report


# Grafisch: ggstatsplot


# Tabellarisch: gtsummary

Korrelation

Standard-Ausgabe

Korrelation <- cor.test(topalbums$raw_usa,
                        topalbums$raw_eur)

Korrelation

    Pearson's product-moment correlation

data:  topalbums$raw_usa and topalbums$raw_eur
t = -3.7065, df = 2998, p-value = 0.0002139
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.10307656 -0.03182962
sample estimates:
        cor 
-0.06753919 

Text-Interpretation

report(Korrelation)

The Pearson's product-moment correlation between topalbums$raw_usa and topalbums$raw_eur is negative, statistically significant, and very small (r = -0.07, 95% CI [-0.10, -0.03], t(2998) = -3.71, p < .001)

Grafisch

ggscatterstats(data = topalbums,
    x = raw_usa, y = raw_eur)

Simple Syntax - informative Ausgabe!

• Streudiagramm: bivariate Verteilungen
• Dichte + Histogramm: “marginal densigram” - Verteilungen
• Anpassungsgerade mit Konfidenzband: Linearer Zusammenhang?
• Korrelationskoeffizient Pearson’s r
• t-Wert, p-Wert, Konfidenzintervall für r, Fallzahl
• Bayes Factor

Zum Üben: Bayes-Faktor verbergen

Tipp: Parameter bf.message; Voreinstellung: TRUE

# ggscatterstats(data = topalbums,
#     x = raw_usa, y = raw_eur,
#     fb.message = xxx)

Regression

• Nimmt im Gegensatz zur Korrelation eine Wirkungsrichtung an
• Unterscheidung zwischen unabhängigen Variablen (Prädiktoren, Regressoren) und abhängiger Variable (Zielgröße, Regressand)

Standard-Ausgabe

reg <- lm(raw_eur ~ raw_usa,
          data = topalbums)

reg

Call:
lm(formula = raw_eur ~ raw_usa, data = topalbums)

Coefficients:
(Intercept)      raw_usa  
    4.51192     -0.07867  

Hier ist das summary so programmiert, dass es eine ausführlichere Ausgabe enthält als das Regressionsobjekt selbst.

summary(reg)

Call:
lm(formula = raw_eur ~ raw_usa, data = topalbums)

Residuals:
   Min     1Q Median     3Q    Max 
-4.403 -2.923 -1.079  1.761 14.900 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.51192    0.11626  38.809  < 2e-16 ***
raw_usa     -0.07867    0.02122  -3.707 0.000214 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.793 on 2998 degrees of freedom
Multiple R-squared:  0.004562,  Adjusted R-squared:  0.00423 
F-statistic: 13.74 on 1 and 2998 DF,  p-value: 0.0002139

Text-Interpretation: report

# Zum Selbst-Ausfüllen ...

Grafische Darstellung siehe Korrelation.

Test auf Normalverteilung

Voraussetzung vieler parametrischer Verfahren wie z. B. t-Test.

Statistisch

shapiro.test(topalbums$final_score)

    Shapiro-Wilk normality test

data:  topalbums$final_score
W = 0.6974, p-value < 2.2e-16

Anregung zum Üben: Andere Variable auswählen, z. B. year oder andere Punktzahl wie raw_usa, oder anderen Datensatz verwenden (z. B. songs2000 mit indicativerevenue)

Grafisch: ggstatsplot

ggstatsplot::gghistostats(
  topalbums, final_score,
  normal.curve = TRUE)

Anregungen zum Üben:

• Andere Effektgröße auswählen. Voreinstellung: Hedge’s g (“g”); Alternative: Cohen’s d (“d”)

# ggstatsplot::gghistostats(
#   topalbums, final_score,
#   normal.curve = TRUE,
#   effsize.type = "xxx")

Ergänzungen

Interaktive Grafiken

Ein “normales” statisches Diagramm …

Anstatt den ggplot-Befehl direkt anzuwenden, kann man das Diagramm einem Objekt zuweisen: p für plot.

Dieses können wir dann später weiterverwenden.

p <- topalbums %>% 
  filter(artist %in% c("Bob Dylan", "The Beatles", "The Rolling Stones", "Bruce Springsteen", "Madonna")) %>% 
  ggplot(aes(x = year, y = final_score, label = name, col = artist)) +
  labs(title = "Alben & Punktzahlen",
       subtitle = "Ausgewählte Bands",
       x = "Jahr", y = "Punktzahl") +
  scale_color_discrete(name = "Künstler / Band") +
  geom_point()

p

# Kurz für: print(p)

Was müssen wir tun, um dieses Diagramm interaktiv zu machen?

plotly::ggplotly(p)

Anregung zum Üben: Andere Bands auswählen; wer mag, kann einen anderen Datensatz verwenden (z. B. songs2000 oder albums2000).

Animationen mit gganimate

Code von Thomas L. Pedersen, siehe Dokumentation zu gganimate.

library(gapminder)
library(ggplot2)
library(gganimate)

ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, colour = country)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = country_colors) +
  scale_size(range = c(2, 12)) +
  scale_x_log10() +
  facet_wrap(~continent) +
  theme_bw() +
  theme(legend.position = 'none') +
  labs(title = 'Year: {frame_time}', x = 'GDP per capita', y = 'life expectancy') +
  transition_time(year)

Regressionstabellen mit texreg::htmlreg

Ein Modell:

Stones_lm <- topalbums %>% 
  filter(artist == "The Rolling Stones") %>% 
  lm(raw_eur ~ raw_usa, data = .)

htmlreg(Stones_lm, doctype = FALSE)

Statistical models

	Model 1
(Intercept)	5.45***
	(1.22)
raw_usa	-0.09
	(0.18)
R2	0.01
Adj. R2	-0.02
Num. obs.	36
***p < 0.001; **p < 0.01; *p < 0.05

htmlreg(Stones_lm, single.row = TRUE, doctype = FALSE)

Statistical models

	Model 1
(Intercept)	5.45 (1.22)***
raw_usa	-0.09 (0.18)
R2	0.01
Adj. R2	-0.02
Num. obs.	36
***p < 0.001; **p < 0.01; *p < 0.05

Vergleich von zwei Modellen:

Beatles_lm <- topalbums %>%
  filter(artist == "The Beatles") %>% 
  lm(raw_eur ~ raw_usa, data = .)

htmlreg(list(Stones_lm, Beatles_lm),
        custom.model.names = c("Stones", "Beatles"),
        doctype = FALSE)

Statistical models

	Stones	Beatles
(Intercept)	5.45***	3.18**
	(1.22)	(1.03)
raw_usa	-0.09	0.18
	(0.18)	(0.09)
R2	0.01	0.12
Adj. R2	-0.02	0.09
Num. obs.	36	32
***p < 0.001; **p < 0.01; *p < 0.05

Das gt-Paket für flexibel gestaltbare Tabellen

Das gt-Paket ermöglicht sehr flexibel gestaltbare Tabellen mit relativ intuitiver Syntax.

topalbums %>% 
  head() %>% 
  gt()

position	artist	name	year	final_score	raw_usa	raw_eng	raw_eur	ueber5
1	The Beatles	Sgt Pepper's Lonely Hearts Club Band	1967	54.415	33.740	13.824	9.284	> 5 P
2	Michael Jackson	Thriller	1982	50.782	29.628	18.143	14.952	> 5 P
3	U2	The Joshua Tree	1987	38.197	21.342	13.411	17.733	> 5 P
4	The Beatles	Abbey Road	1969	37.233	24.587	9.409	9.953	> 5 P
5	Fleetwood Mac	Rumours	1977	36.616	30.057	15.097	7.300	> 5 P
6	Pink Floyd	Dark Side Of The Moon	1973	35.338	26.903	13.621	9.717	> 5 P

Nehmen wir an, final_score wäre in US-Dollar - wie können wir die Zahlen einfach formatieren?

Dazu Tabellenüberschrift und Unter-Überschrift; Zellen mit raw_usa > 24 hervorgehoben.

topalbums %>%
  head() %>% 
  gt() %>% 
  fmt_currency(columns = "final_score",
               currency = "USD") %>% 
  tab_header(title = "Top Alben",
             subtitle = glue::glue("Von {min(topalbums$year)} bis {max(topalbums$year)}")) %>% 
  tab_style(
    style = list(
      cell_fill(color = "lightblue"),
      cell_text(style = "italic")
    ),
    locations = cells_body(
      columns = raw_usa,
      rows = raw_usa > 24
    )
  )

Top Alben
Von 1940 bis 2018
position	artist	name	year	final_score	raw_usa	raw_eng	raw_eur	ueber5
1	The Beatles	Sgt Pepper's Lonely Hearts Club Band	1967	$54.41	33.740	13.824	9.284	> 5 P
2	Michael Jackson	Thriller	1982	$50.78	29.628	18.143	14.952	> 5 P
3	U2	The Joshua Tree	1987	$38.20	21.342	13.411	17.733	> 5 P
4	The Beatles	Abbey Road	1969	$37.23	24.587	9.409	9.953	> 5 P
5	Fleetwood Mac	Rumours	1977	$36.62	30.057	15.097	7.300	> 5 P
6	Pink Floyd	Dark Side Of The Moon	1973	$35.34	26.903	13.621	9.717	> 5 P